/* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using System;
using System.Collections.Generic;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra
{
    public partial class SparseMatrix
    {
        // UPGRADE NOTE: Could change the current storage scheme from
        // Compressed Row Storage (currently implemented) to a better one
        // The literature has several options. Suggestions are:
        // 1) a special row storage scheme which uses individual arrays for each 
        //    row. There would be one object array for the row 'list'. Each item 
        //    in the row array would contain either null (no values in that row)
        // 2) An array of structs (the struct contains the value and the column index) 
        //    which acts as a column array. Only non zero values are stored
        //    This way rowwise access may be a bit quicker and as the columns 
        //    will in general be smaller than a full column they might also
        //    fit in the CPU cache better.
        // 3) Use a compressed storage. Several papers describe the use of
        //    a compressed storage which will reduce memory bandwidth.
        //
        // Both 1) and 2) should result in smaller arrays being used which could mean
        // we don't allocate on the Large Object Heap
        //
        // Before converting any to any of these storage systems we'll have to write 
        // some tests that can tell us which is faster. Suggested tests are:
        // 1) Get value
        // 2) Set value
        // 3) Copy matrix
        // 4) Copy vector
        // 5) Arithmetic operations (+, -, *)
        // 6) Solver application

        // Can also use sorted lists for the sub columns.

        /// <summary>
        /// The array containing the actual non-zero values. Stored in a compressed
        /// row fashion.
        /// </summary>
        internal double[] mValues = new double[0];

        /// <summary>
        /// An array containing the column indices of the non-zero values.
        /// </summary>
        internal int[] mColumnIndices = new int[0];

        /// <summary>
        /// The array containing the row indices of the existing rows. 
        /// </summary>
        internal int[] mRowIndices = new int[0];

        /// <summary>
        /// The total number of non-zero values in the matrix.
        /// </summary>
        /// <remarks>
        /// This variable keeps track of the number of non-zero values that are stored
        /// in the <b>mValues</b> and <b>mColumnIndices</b> array's. Because of this
        /// <b>mValueCount</b> also indicates the 'length' of the array's. Note that the
        /// actual arrays may be longer (but never shorter) than indicated by
        /// <b>mValueCount</b>. If the actual arrays are longer than indicated reading
        /// beyond the number of values (i.e. <b>mValueCount</b>) will produce garbage
        /// data.
        /// </remarks>
        internal int mValueCount;

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <param name="rows">The number of rows in the matrix.</param>
        /// <param name="columns">The number of columns in the matrix.</param>
        public SparseMatrix(int rows, int columns) : base(rows, columns)
        {
            mRowIndices = new int[rows];
        }

        #region Locate matrix item in the arrays

        // UPGRADE NOTE: Could potentially store last accessed item (the index numbers
        // because most of the time matrices are accessed sequentially which means
        // that we can figure out what item to pick without having to actually
        // search for it.

        /// <summary>
        /// Finds either the item indicated by the given indices (returns <b>true</b>) or the 
        /// insert position (returns <b>false</b>) in the values array. 
        /// </summary>
        /// <param name="rowIndex"></param>
        /// <param name="columnIndex"></param>
        /// <param name="itemIndex"></param>
        /// <returns></returns>
        private bool FindItem(int rowIndex, int columnIndex, out int itemIndex)
        {
            if ((rowIndex < 0) || (rowIndex >= Rows) || (columnIndex < 0)
                || (columnIndex >= Columns))
            {
                itemIndex = -1;
                return false;
            }

            // First find the start and end indices for the row
            int startIndex;
            int endIndex;
            FindRowStartAndEnd(rowIndex, out startIndex, out endIndex);

            // Check if startIndex and endIndex are equal. If so then
            // the item does not exist
            if (startIndex == endIndex)
            {
                itemIndex = endIndex;
                return false;
            }

            // Now subtract 1 from the endIndex because the current endIndex
            // indicates were the next row starts. 
            endIndex -= 1;

            // If the startIndex is only one bigger than the endIndex
            // do a separate check because the binary search loop will not 
            // be entered.
            // Check if one of the boundaries is the correct index
            if (startIndex >= endIndex - 1)
            {
                return FindItemForSizeTwoArray(startIndex, endIndex, columnIndex, out itemIndex);
            }

            itemIndex = (startIndex + endIndex) / 2;
            // Start the binary search
            while (startIndex < endIndex - 1)
            {
                if (mColumnIndices[itemIndex] == columnIndex)
                {
                    // Found the value we are looking for so we can exit
                    // here
                    return true;
                }

                // We're still here so we didn't find it. Now see in
                // what direction to search
                if (mColumnIndices[itemIndex] > columnIndex)
                {
                    // value is to big that means that we need to move 
                    // back in the array
                    endIndex = itemIndex;
                }
                else
                {
                    // So the value must be to small
                    // Move forward in the array
                    startIndex = itemIndex;
                }
                itemIndex = (startIndex + endIndex) / 2;

                // Do a quick check to see if the end of the search was
                // reached. This occurs when itemIndex is equal to one
                // of the two boundaries. In this case this should always
                // be the lower end because integer division truncates instead
                // of rounding
                if (startIndex >= endIndex - 1)
                {
                    return FindItemForSizeTwoArray(startIndex, endIndex, columnIndex, out itemIndex);
                }
            }
            return false;
        }

        /// <summary>
        /// Locates the start and end of the row with the specified index. 
        /// </summary>
        /// <param name="rowIndex">The index of the row that must be located.</param>
        /// <param name="startIndex">The index of the first element in the <see cref="mValues"/> array that belongs
        /// to the given row.</param>
        /// <param name="endIndex">The index of the first element in the <see cref="mValues"/> array that belongs
        /// to the given row + 1.</param>
        private void FindRowStartAndEnd(int rowIndex, out int startIndex, out int endIndex)
        {
            // First find the start and end indices for the row
            startIndex = mRowIndices[rowIndex];
            if (rowIndex < mRowIndices.Length - 1)
            {
                endIndex = mRowIndices[rowIndex + 1];
            }
            else
            {
                // We are looking for a value in the last row so
                // the end of the search is reached as we get to the
                // last valid item in the column array
                endIndex = mValueCount;
            }
        }

        private bool FindItemForSizeTwoArray(int startIndex, int endIndex, int columnIndex, out int itemIndex)
        {
            if ((startIndex > -1) && (startIndex < mColumnIndices.Length)
                &&
                (mColumnIndices[startIndex] == columnIndex))
            {
                itemIndex = startIndex;
                return true;
            }

            if ((endIndex > -1) && (endIndex < mColumnIndices.Length)
                &&
                (mColumnIndices[endIndex] == columnIndex))
            {
                itemIndex = endIndex;
                return true;
            }

            // neither one of the boundaries is actually correct
            // Check where eventual new values should be placed
            if (mColumnIndices[startIndex] > columnIndex)
            {
                itemIndex = startIndex;
            }
            else
            {
                if (mColumnIndices[endIndex] > columnIndex)
                {
                    itemIndex = endIndex;
                }
                else
                {
                    itemIndex = endIndex + 1;
                }
            }
            return false;
        }

        #endregion

        #region Resize the data storage arrays

        private void GrowStorage()
        {
            // If we already are at the maximum size then exit now
            if (mValues.Length == (Rows * Columns))
            {
                return;
            }

            // Determine the increase in size. We will not grow beyond
            // the size of the matrix
            int size = System.Math.Min(mValues.Length + GrowthSize(), Rows * Columns);

            // Grow the valueArray by X%
            // Grow the columnIndex array by X%
            // Both arrays should always be equal in length!!!!!
            double[] newValues = new double[size];
            int[] newIndices = new int[size];
            
            Buffer.BlockCopy(mValues, 0, newValues, 0, mValues.Length * Constants.SizeOfDouble);
            Buffer.BlockCopy(mColumnIndices, 0, newIndices, 0, mColumnIndices.Length * Constants.SizeOfInt);

            mValues = newValues;
            mColumnIndices = newIndices;
        }

        private int GrowthSize()
        {
            int delta;
            if (mValues.Length > 1024)
            {
                delta = mValues.Length / 4;
            }
            else
            {
                if (mValues.Length > 256)
                {
                    delta = 512;
                }
                else
                {
                    if (mValues.Length > 64)
                    {
                        delta = 128;
                    }
                    else
                    {
                        delta = 32;
                    }
                }
            }
            return delta;
        }

        #endregion

        #region Internal add item to the arrays

        private void InsertValue(int rowIndex, int columnIndex, double value)
        {
            // Check if the value is not zero.
            if (!Precision.EqualsWithinDecimalPlaces(0.0, value, Constants.DecimalsForComparison))
            {
                // We assume that the item is not present in the array
                // First find the spot where the value must be inserted
                int index;
                FindItem(rowIndex, columnIndex, out index);

                // Check if the storage needs to be increased
                if (mValueCount == mValues.Length)
                {
                    // Value array is completely full so we increase the size
                    // This also grows the columIndex array
                    GrowStorage();
                }
                // Move all values (with an position larger than index) in the value array 
                // to the next position
                // move all values (with an position larger than index) in the columIndices 
                // array to the next position
                for (int i = mValueCount - 1; i > index - 1; i--)
                {
                    mValues[i + 1] = mValues[i];
                    mColumnIndices[i + 1] = mColumnIndices[i];
                }

                // Add the value and the column index
                mValues[index] = value;
                mColumnIndices[index] = columnIndex;

                // increase the number of non-zero numbers by one
                mValueCount += 1;

                // add 1 to all the row indices for rows bigger than rowIndex
                // so that they point to the correct part of the value array 
                // again.
                for (int i = rowIndex + 1; i < mRowIndices.Length; i++)
                {
                    mRowIndices[i] += 1;
                }
            }
        }

        #endregion

        #region Internal remove item from the arrays

        private void DeleteItem(int rowIndex, int columnIndex)
        {
            // We have an item that needs to be removed from the
            // array. First check if it is actually stored in the array
            int index;
            if (FindItem(rowIndex, columnIndex, out index))
            {
                // Move all values (with an position larger than index) in the value array 
                // to the previous position
                // move all values (with an position larger than index) in the columIndices 
                // array to the previous position
                for (int i = index + 1; i < mValueCount; i++)
                {
                    mValues[i - 1] = mValues[i];
                    mColumnIndices[i - 1] = mColumnIndices[i];
                }

                for (int i = rowIndex + 1; i < mRowIndices.Length; i++)
                {
                    mRowIndices[i] -= 1;
                }
                mValueCount -= 1;
            }
        }

        #endregion

        /// <summary>
        /// Retrieves the value at the given indices. <strong>This method is not ranged checked.</strong> If the
        /// indices are outside the dimensions of the matrix, it may return an incorrect value or throw
        /// an <see cref="IndexOutOfRangeException"/>.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <returns>The value at the given indices.</returns>
        /// <exception cref="IndexOutOfRangeException">Depending on the implementation, an <see cref="IndexOutOfRangeException"/>
        /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
        protected internal override double ValueAt(int row, int column)
        {
            int index;
            if (FindItem(row, column, out index))
            {
                return mValues[index];
            }
            return 0.0;
        }

        /// <summary>
        /// Set the value at the given indices to the given value. <strong>This method is not ranged checked.</strong> If the
        /// indices are outside the dimensions of the matrix, it may set an incorrect element or throw
        /// an <see cref="IndexOutOfRangeException"/>.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <param name="value">The value to set.</param>
        /// <exception cref="IndexOutOfRangeException">Depending on the implementation, an <see cref="IndexOutOfRangeException"/>
        /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
        protected internal override void ValueAt(int row, int column, double value)
        {
            int index;
            if (FindItem(row, column, out index))
            {
                // Check if the value is zero.
                if (Precision.EqualsWithinDecimalPlaces(0.0, value, Constants.DecimalsForComparison))
                {
                    DeleteItem(row, column);
                }
                else
                {
                    mValues[index] = value;
                }
            }
            else
            {
                InsertValue(row, column, value);
            }
        }

        /// <summary>
        /// Returns the number of non zero elements in the matrix.
        /// </summary>
        /// <value>The number of non zero elements.</value>
        public int NonZeros
        {
            get
            {
                return mValueCount;
            }
        }

        /// <summary>Calculates the L1 norm.</summary>
        /// <returns>The L1 norm of the matrix.</returns>
        public override double L1Norm()
        {
            // Allocate an array which stores the intermediate results
            var sums = new double[Columns];

            // The L1 norm is the largest column sum of the matrix
            for (int i = 0; i < Rows; i++)
            {
                int startIndex;
                int endIndex;
                FindRowStartAndEnd(i, out startIndex, out endIndex);

                for (int j = startIndex; j < endIndex; j++)
                {
                    sums[mColumnIndices[j]] += System.Math.Abs(mValues[j]);
                }
            }

            double norm = double.NegativeInfinity;
            for (int i = 0; i < sums.Length; i++)
            {
                norm = System.Math.Max(norm, sums[i]);
            }

            return norm;
        }

        /// <summary>Calculates the L2 norm.</summary>
        /// <returns>The L2 norm of the matrix.</returns>   
        /// <remarks>For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition. 
        /// In a later release, it will be replaced with a sparse implementation.</remarks>
        public override double L2Norm()
        {
            // TODO: Implement a proper sparse version of the L2Norm
            return base.L2Norm();
        }

        /// <summary>Calculates the Frobenius norm of this matrix.</summary>
        /// <returns>The Frobenius norm of this matrix.</returns>
        public override double FrobeniusNorm()
        {
            // norm = sqrt(the sum of the squares of all items)
            double norm = 0;
            for (int i = 0; i < mValueCount; i++)
            {
                norm += mValues[i] * mValues[i];
            }

            return System.Math.Sqrt(norm);
        }

        /// <summary>Calculates the infinity norm of this matrix.</summary>
        /// <returns>The infinity norm of this matrix.</returns>      
        public override double InfinityNorm()
        {
            // The infinity norm is the largest row sum of the matrix
            // This one is simple because all is sorted by rows
            double norm = double.NegativeInfinity;
            for (int i = 0; i < Rows; i++)
            {
                int startIndex;
                int endIndex;
                FindRowStartAndEnd(i, out startIndex, out endIndex);

                double rowSum = 0;
                for (int j = startIndex; j < endIndex; j++)
                {
                    rowSum += System.Math.Abs(mValues[j]);
                }
                norm = System.Math.Max(norm, rowSum);
            }
            return norm;
        }

        /// <summary>Calculates the condition number of this matrix.</summary>
        /// <returns>The condition number of the matrix.</returns>
        /// <remarks>The condition number is calculated using singular value decomposition.</remarks>
        /// <exception cref="ConvergenceFailedException">If the singular value decomposition
        /// cannot converge to a solution.</exception>
        public override double ConditionNumber()
        {
            // @TODO: Implement a proper sparse version of the ConditionNumber
            return base.ConditionNumber();
        }

        /// <summary>
        /// Sets all values to zero.
        /// </summary>
        public override void Clear()
        {
            mValueCount = 0;
            Array.Clear(mRowIndices, 0, mRowIndices.Length);
        }

        /// <summary>
        /// Copies a column into to the given vector.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <param name="result">The vector to copy the column into.</param>
        /// <exception cref="ArgumentNullException">If the result matrix is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of columns.</exception>        
        /// <exception cref="NotConformableException">If the result vector is of the incorrect size.</exception>
        public override void GetColumn(int index, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (index >= Columns || index < 0)
            {
                throw new ArgumentOutOfRangeException("index");
            }

            if (result.Count != Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the vector
            result.Clear();

            for (int i = 0; i < Rows; i++)
            {
                int itemIndex;
                if (!FindItem(i, index, out itemIndex))
                {
                    continue;
                }

                result[i] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies the requested column elements into the given vector.
        /// </summary>
        /// <param name="columnIndex">The column to copy elements from.</param>
        /// <param name="rowIndex">The row to start copying from.</param>
        /// <param name="length">The number of elements to copy.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the column into.</param>
        /// <exception cref="ArgumentNullException">If the result <see cref="Vector"/> is <see langword="null" />.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is negative,
        /// or greater than or equal to the number of columns.</exception>        
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is negative,
        /// or greater than or equal to the number of rows.</exception>        
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> + <paramref name="length"/>  
        /// is greater than or equal to the number of rows.</exception>
        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
        /// <exception cref="NotConformableException">If <strong>result.Count &lt; length</strong>.</exception>
        public override void GetColumn(int columnIndex, int rowIndex, int length, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (columnIndex >= Columns || columnIndex < 0)
            {
                throw new ArgumentOutOfRangeException("columnIndex");
            }
            if (rowIndex >= Rows || rowIndex < 0)
            {
                throw new ArgumentOutOfRangeException("rowIndex");
            }

            if (rowIndex + length > Rows)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            if (length < 1)
            {
                throw new ArgumentException(Resources.NotPositive, "length");
            }

            if (result.Count < length)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result vector
            result.Clear();

            for (int i = rowIndex, k = 0; i < rowIndex + length; i++, k++)
            {
                int itemIndex;
                if (!FindItem(i, columnIndex, out itemIndex))
                {
                    continue;
                }

                result[k] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies a row into to the given <see cref="Vector"/>.
        /// </summary>
        /// <param name="index">The row to copy.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the row into.</param>
        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of rows.</exception>
        /// <exception cref="NotConformableException">If <b>this.Columns != result.Count</b>.</exception>
        public override void GetRow(int index, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (index >= Rows
                || index < 0)
            {
                throw new ArgumentOutOfRangeException("index");
            }

            if (result.Count != Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result vector
            result.Clear();

            // The data is stored in compressed row storage so all we have to do is find the
            // begin and end of the row and we can quickly iterate over that.
            int begin;
            int end;
            FindRowStartAndEnd(index, out begin, out end);

            for (int i = begin; i < end; i++)
            {
                int insertPosition = mColumnIndices[i];
                result[insertPosition] = mValues[i];
            }
        }

        /// <summary>
        /// Copies the requested row elements into a new <see cref="Vector"/>.
        /// </summary>
        /// <param name="rowIndex">The row to copy elements from.</param>
        /// <param name="columnIndex">The column to start copying from.</param>
        /// <param name="length">The number of elements to copy.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the column into.</param>
        /// <exception cref="ArgumentNullException">If the result <see cref="Vector"/> is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is negative,
        /// or greater than or equal to the number of columns.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is negative,
        /// or greater than or equal to the number of rows.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> + <paramref name="length"/>
        /// is greater than or equal to the number of rows.</exception>
        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
        /// <exception cref="NotConformableException">If <strong>result.Count &lt; length</strong>.</exception>
        public override void GetRow(int rowIndex, int columnIndex, int length, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (rowIndex >= Rows || rowIndex < 0)
            {
                throw new ArgumentOutOfRangeException("rowIndex");
            }
            if (columnIndex >= Columns || columnIndex < 0)
            {
                throw new ArgumentOutOfRangeException("columnIndex");
            }

            if (columnIndex + length > Columns)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            if (length < 1)
            {
                throw new ArgumentException(Resources.NotPositive, "length");
            }

            if (result.Count < length)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result vector
            result.Clear();

            for (int i = columnIndex, k = 0; i < columnIndex + length; i++, k++)
            {
                int itemIndex;
                if (!FindItem(rowIndex, i, out itemIndex))
                {
                    continue;
                }

                result[k] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies the elements of this matrix to the given matrix.
        /// </summary>
        /// <param name="target">The matrix to copy values into.</param>
        /// <exception cref="ArgumentNullException">If target is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If target is not the same size as this matrix.</exception>
        public override void CopyTo(Matrix target)
        {
            if (target == null)
            {
                throw new ArgumentNullException("target");
            }

            if ((target.Rows != Rows) || (target.Columns != Columns))
            {
                throw new NotConformableException("target", Resources.ParameterNotConformable);
            }

            if( ReferenceEquals(this, target))
            {
                return;
            }

            SparseMatrix matrix = target as SparseMatrix;
            if (matrix != null)
            {
                // The matrix is a sparse matrix. Do a quick block copy of the
                // data.
                CopyToSparseMatrix(matrix);
            }
            else
            {
                // Clear the target matrix
                target.Clear();

                // Do a normal copy
                // iterate over the rows
                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    // Get the begin / end index for the current row
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    // Get the values for the current row
                    if (begin == end)
                    {
                        // Begin and end are equal. There are no values in the row
                        // Move to the next row
                        continue;
                    }

                    for (int j = begin; j < end; j++)
                    {
                        // We could use the ValueAt method.
                        // target.ValueAt(i, mColumnIndices[j], mValues[j]);

                        target[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        private void CopyToSparseMatrix(SparseMatrix matrix)
        {
            // copy the row indices
            Buffer.BlockCopy(mRowIndices, 0, matrix.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);

            // copy the column indices. Only copy the actual data items. Note that the array
            // may be bigger than the number of items we copy.
            matrix.mColumnIndices = new int[mValueCount];
            Buffer.BlockCopy(mColumnIndices, 0, matrix.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);

            // copy the data. Only copy the actual data items. Note that the array
            // may be bigger than the number of items we copy.
            matrix.mValueCount = mValueCount;
            matrix.mValues = new double[mValueCount];
            Buffer.BlockCopy(mValues, 0, matrix.mValues, 0, mValueCount * Constants.SizeOfDouble);
        }

        /// <summary>Computes the determinant of this matrix.</summary>
        /// <returns>The determinant of this matrix.</returns>
        /// <exception cref="MatrixNotSquareException">If the matrix is not square.</exception>
        /// <remarks>The determinant is calculated using LU decomposition.</remarks>
        public override double Determinant()
        {
            // TODO: Implement proper sparse version of Determinant
            return base.Determinant();
        }
        
        /// <summary>
        /// Returns the transpose of this matrix.
        /// </summary>
        /// <returns>The transpose of this matrix.</returns>
        public override Matrix Transpose()
        {
            Matrix result = CreateMatrix(Columns, Rows);

            // Do an 'inverse' CopyTo
            // iterate over the rows
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                // Get the begin / end index for the current row
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                // Get the values for the current row
                if (begin == end)
                {
                    // Begin and end are equal. There are no values in the row
                    // Move to the next row
                    continue;
                }

                for (int j = begin; j < end; j++)
                {
                    // We could use the ValueAt method.
                    // target.ValueAt(i, mColumnIndices[j], mValues[j]);

                    result[mColumnIndices[j], i] = mValues[j];
                }
            }

            return result;
        }

        /// <summary>
        /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
        /// </summary>
        /// <returns>The inverse of this matrix.</returns>
        /// <exception cref="MatrixNotSquareException">If this matrix is not square.</exception>
        public override Matrix Inverse()
        {
            // TODO: implement a proper sparse version of the Inverse method
            return base.Inverse();
        }

        /// <summary>
        /// Returns a new matrix containing the lower triangle of this matrix.
        /// </summary>
        /// <returns>The lower triangle of this matrix.</returns>        
        public override Matrix GetLowerTriangle()
        {
            Matrix result = CreateMatrix(Rows, Columns);

            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i >= mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Puts the lower triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the same as this matrix.</exception>
        public override void GetLowerTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                GetLowerTriangle(tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result matrix
                result.Clear();

                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    for (int j = begin; j < end; j++)
                    {
                        if (i >= mColumnIndices[j])
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
        /// does not contain the diagonal elements of this matrix.
        /// </summary>
        /// <returns>The lower triangle of this matrix.</returns>
        public override Matrix GetStrictlyLowerTriangle()
        {
            Matrix result = CreateMatrix(Rows, Columns);
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i > mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
            return result;
        }

        /// <summary>
        /// Puts the strictly lower triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the same as this matrix.</exception>
        public override void GetStrictlyLowerTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                GetStrictlyLowerTriangle(tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result matrix
                result.Clear();

                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    for (int j = begin; j < end; j++)
                    {
                        if (i > mColumnIndices[j])
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Returns a new matrix containing the upper triangle of this matrix.
        /// </summary>
        /// <returns>The upper triangle of this matrix.</returns>   
        public override Matrix GetUpperTriangle()
        {
            Matrix result = CreateMatrix(Rows, Columns);
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i <= mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
            return result;
        }

        /// <summary>
        /// Puts the upper triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the same as this matrix.</exception>
        public override void GetUpperTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if ((result.Rows != Rows) || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                GetUpperTriangle(tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result matrix
                result.Clear();

                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    for (int j = begin; j < end; j++)
                    {
                        if (i <= mColumnIndices[j])
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Returns a new matrix containing the upper triangle of this matrix. The new matrix
        /// does not contain the diagonal elements of this matrix.
        /// </summary>
        /// <returns>The upper triangle of this matrix.</returns>
        public override Matrix GetStrictlyUpperTriangle()
        {
            Matrix result = CreateMatrix(Rows, Columns);
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i < mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
            return result;
        }

        /// <summary>
        /// Puts the strictly upper triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the same as this matrix.</exception>
        public override void GetStrictlyUpperTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                GetStrictlyUpperTriangle(tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result matrix
                result.Clear();

                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    for (int j = begin; j < end; j++)
                    {
                        if (i < mColumnIndices[j])
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Creates a matrix that contains the values from the requested sub-matrix.
        /// </summary>
        /// <param name="rowIndex">The row to start copying from.</param>
        /// <param name="rowLength">The number of rows to copy. Must be positive.</param>
        /// <param name="columnIndex">The column to start copying from.</param>
        /// <param name="columnLength">The number of columns to copy. Must be positive.</param>
        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
        /// negative, or greater than or equal to the number of rows.</item>
        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
        /// of columns.</item>
        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
        /// <exception cref="ArgumentException">If <paramref name="rowLength"/> or <paramref name="columnLength"/>
        /// is not positive.</exception>
        public override Matrix GetSubMatrix(int rowIndex, int rowLength, int columnIndex, int columnLength)
        {
            if (rowIndex >= Rows || rowIndex < 0)
            {
                throw new ArgumentOutOfRangeException("rowIndex");
            }

            if (columnIndex >= Columns || columnIndex < 0)
            {
                throw new ArgumentOutOfRangeException("columnIndex");
            }

            if (rowLength < 1)
            {
                throw new ArgumentException(Resources.NotPositive, "rowLength");
            }
            if (columnLength < 1)
            {
                throw new ArgumentException(Resources.NotPositive, "columnLength");
            }
            int colMax = columnIndex + columnLength;
            int rowMax = rowIndex + rowLength;

            if (rowMax > Rows)
            {
                throw new ArgumentOutOfRangeException("rowLength");
            }

            if (colMax > Columns)
            {
                throw new ArgumentOutOfRangeException("columnLength");
            }

            Matrix result = CreateMatrix(rowLength, columnLength);

            for (int i = rowIndex, row = 0; i < rowIndex + rowLength; i++, row++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    // check if the column index is in the range
                    if ((mColumnIndices[j] >= columnIndex) && (mColumnIndices[j] < columnIndex + columnLength))
                    {
                        int column = mColumnIndices[j] - columnIndex;
                        result[row, column] = mValues[j];
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Returns this matrix as a multidimensional array.
        /// </summary>
        /// <returns>a multidimensional containing the values of this matrix.</returns>        
        public override double[,] ToArray()
        {
            double[,] result = new double[Rows, Columns];
            for (int rowIndex = 0; rowIndex < mRowIndices.Length; rowIndex++)
            {
                int columnStart;
                int columnEnd;
                FindRowStartAndEnd(rowIndex, out columnStart, out columnEnd);

                for (int columnIndex = columnStart; columnIndex < columnEnd; columnIndex++)
                {
                    result[rowIndex, mColumnIndices[columnIndex]] = mValues[columnIndex];
                }
            }
            return result;
        }

        /// <summary>
        /// Adds another matrix to this matrix.
        /// </summary>
        /// <param name="other">The matrix to add to this matrix.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Add(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if ((other.Rows != Rows) || (other.Columns != Columns))
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, other))
            {
                Multiply(2);
                return;
            }

            SparseMatrix matrix = other as SparseMatrix;
            if (matrix != null)
            {
                AddSparseMatrixToSelf(matrix);
            }
            else
            {
                for (int i = 0; i < other.Rows; i++)
                {
                    for (int j = 0; j < other.Columns; j++)
                    {
                        int index;
                        if (FindItem(i, j, out index))
                        {
                            // check if a is the negative of mValues[index].
                            // If so then we just delete the value at index
                            if (Precision.EqualsWithinDecimalPlaces(0.0, mValues[index] + other[i, j], Constants.DecimalsForComparison))
                            {
                                DeleteItem(i, j);
                            }
                            else
                            {
                                mValues[index] += other[i, j];
                            }
                        }
                        else
                        {
                            InsertValue(i, j, other[i, j]);
                        }
                    }
                }
            }
        }

        private void AddSparseMatrixToSelf(SparseMatrix matrix)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                int begin;
                int end;
                matrix.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    int index;
                    if (FindItem(i, matrix.mColumnIndices[j], out index))
                    {
                        // check if a is the negative of mValues[index].
                        // If so then we just delete the value at index
                        if (Precision.EqualsWithinDecimalPlaces(0.0, mValues[index] + matrix.mValues[j], Constants.DecimalsForComparison))
                        {
                            DeleteItem(i, matrix.mColumnIndices[j]);
                        }
                        else
                        {
                            mValues[index] += matrix.mValues[j];
                        }
                    }
                    else
                    {
                        InsertValue(i, matrix.mColumnIndices[j], matrix.mValues[j]);
                    }
                }
            }
        }

        /// <summary>
        /// Negates each element of this matrix.
        /// </summary>        
        public override void Negate()
        {
            for (int i = 0; i < mValueCount; i++)
            {
                mValues[i] = -mValues[i];
            }
        }

        /// <summary>
        /// Subtracts another matrix from this matrix.
        /// </summary>
        /// <param name="other">The matrix to subtract.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Subtract(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if ((other.Rows != Rows) || (other.Columns != Columns))
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, other))
            {
                Clear();
                return;
            }

            SparseMatrix matrix = other as SparseMatrix;
            if (matrix != null)
            {
                SubstractSparseMatrixFromSelf(matrix);
            }
            else
            {
                for (int i = 0; i < other.Rows; i++)
                {
                    for (int j = 0; j < other.Columns; j++)
                    {
                        int index;
                        if (FindItem(i, j, out index))
                        {
                            // check if a is the positive of mValues[index].
                            // If so then we just delete the value at index
                            if (Precision.EqualsWithinDecimalPlaces(0.0, mValues[index] - other[i, j], Constants.DecimalsForComparison))
                            {
                                DeleteItem(i, j);
                            }
                            else
                            {
                                mValues[index] -= other[i, j];
                            }
                        }
                        else
                        {
                            InsertValue(i, j, -other[i, j]);
                        }
                    }
                }
            }
        }

        private void SubstractSparseMatrixFromSelf(SparseMatrix matrix)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                int begin;
                int end;
                matrix.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    int index;
                    if (FindItem(i, matrix.mColumnIndices[j], out index))
                    {
                        // check if a is the positive of mValues[index].
                        // If so then we just delete the value at index
                        if (!Precision.EqualsWithinDecimalPlaces(0.0, mValues[index] - matrix.mValues[j], Constants.DecimalsForComparison))
                        {
                            DeleteItem(i, matrix.mColumnIndices[j]);
                        }
                        else
                        {
                            mValues[index] -= matrix.mValues[j];
                        }
                    }
                    else
                    {
                        InsertValue(i, matrix.mColumnIndices[j], -matrix.mValues[j]);
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies each element of this matrix with a scalar overwriting the values of this matrix.
        /// </summary>
        /// <param name="scalar">the scalar to multiply with.</param>
        public override void Multiply(double scalar)
        {
            // if double = 0.0 then just reset the
            // value count to zero
            // Check if the value is not zero.
            if (Precision.EqualsWithinDecimalPlaces(0.0, scalar, Constants.DecimalsForComparison))
            {
                Clear();
                return;
            }

            if (Precision.EqualsWithinDecimalPlaces(1.0, scalar, Constants.DecimalsForComparison))
            {
                return;
            }

            for (int i = 0; i < mValueCount; i++)
            {
                mValues[i] *= scalar;
            }
        }

        /// <summary>
        /// Pointwise multiplies this matrix with another matrix.
        /// </summary>
        /// <param name="other">The matrix to pointwise multiply with this one.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception> 
        /// <exception cref="NotConformableException">If this matrix and <paramref name="other"/> are not the same size.</exception>
        public override Matrix PointwiseMultiply(Matrix other)
        {
            // @TODO: implement proper sparse version of PointwiseMultiply(Matrix)
            return base.PointwiseMultiply(other);
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
        public override void Multiply(Matrix other, Matrix result)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Columns != other.Rows)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            if ((result.Rows != Rows) || (result.Columns != other.Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }


            if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                Multiply(other, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result matrix
                result.Clear();

                // The data is stored in row wise fashion.
                // So compute the values in row wise fashion
                DenseVector columnVector = new DenseVector(other.Rows);

                for (int rowIndex = 0; rowIndex < Rows; rowIndex++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(rowIndex, out begin, out end);
                    for (int columIndex = 0; columIndex < other.Columns; columIndex++)
                    {
                        // clear the vector so that we start with an empty vector every time.
                        // Currently the SparseMatrix.Column doesn't clear the vector.
                        columnVector.Clear();
                        other.GetColumn(columIndex, columnVector);

                        double sum = 0;
                        for (int k = begin; k < end; k++)
                        {
                            sum += mValues[k] * columnVector[mColumnIndices[k]];
                        }
                        result[rowIndex, columIndex] = sum;
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with a vector and places the results into the result matrix.
        /// </summary>
        /// <param name="rightSide">The vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If rightSide is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If <strong>result.Count != this.Rows</strong>.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Columns != rightSide.Count</strong>.</exception>
        public override void Multiply(Vector rightSide, Vector result)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (rightSide.Count != Columns)
            {
                throw new NotConformableException("rightSide", Resources.ParameterNotConformable);
            }

            if (result.Count != Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(rightSide, result))
            {
                Vector tmp = result.CreateVector(result.Count);
                Multiply(rightSide, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result vector
                result.Clear();

                // For each row in the matrix do:
                // 1) find out where the row starts and ends
                // 2) for each column number find the appropriate value in a
                // 3) Multiply a and the value in the matrix
                for (int i = 0; i < Rows; i++)
                {
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    result[i] = 0;
                    for (int j = begin; j < end; j++)
                    {
                        result[i] += mValues[j]*rightSide[mColumnIndices[j]];
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with a vector and places the results into the result matrix. Note the operation is vector * matrix.
        /// </summary>
        /// <param name="leftSide">The vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If leftSide is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <b>result.Count != this.Columns</b>.</exception>
        /// <exception cref="NotConformableException">If <b>this.Rows != leftSide.Count</b>.</exception>
        public override void LeftMultiply(Vector leftSide, Vector result)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (Rows != leftSide.Count)
            {
                throw new NotConformableException("leftSide");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Resources.NullParameterException);
            }

            if (Columns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(leftSide, result))
            {
                Vector tmp = result.CreateVector(result.Count);
                LeftMultiply(leftSide, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                // Clear the result vector
                result.Clear();

                foreach (KeyValuePair<int, Vector> column in GetColumnEnumerator())
                {
                    result[column.Key] = leftSide.DotProduct(column.Value);
                }
            }
        }

        /// <summary>
        /// Stacks this matrix on top of the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="lower">The matrix to stack this matrix upon.</param>
        /// <param name="result">The combined matrix.</param>
        /// <exception cref="ArgumentNullException">If lower is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not (this.Rows + lower.rows) x this.Columns.</exception>
        public override void Stack(Matrix lower, Matrix result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (lower.Columns != Columns)
            {
                throw new NotConformableException("lower", Resources.ParameterNotConformable);
            }

            if ((result.Rows != (Rows + lower.Rows)) || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = lower as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                StackSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    StackSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < lower.Rows; i++)
                    {
                        for (int j = 0; j < lower.Columns; j++)
                        {
                            result[i + Rows, j] = lower[i, j];
                        }
                    }
                }
            }
        }

        private void StackSparseMatricesIntoNewSparseMatrix(SparseMatrix lower, SparseMatrix result)
        {
            // Have all sparse matrices. Just concatenate the values array's
            result.mValues = new double[mValueCount + lower.mValueCount];
            result.mValueCount = mValueCount + lower.mValueCount;
            result.mColumnIndices = new int[mValueCount + lower.mValueCount];


            // Copy the row data
            Buffer.BlockCopy(mRowIndices, 0, result.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mRowIndices.Length; i++)
            {
                result.mRowIndices[i + mRowIndices.Length] = mValueCount + lower.mRowIndices[i];
            }

            // Copy the column data
            Buffer.BlockCopy(mColumnIndices, 0, result.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);
            Buffer.BlockCopy(lower.mColumnIndices, 0, result.mColumnIndices, mValueCount * Constants.SizeOfInt,
                             lower.mValueCount * Constants.SizeOfInt);

            // Copy the values
            Buffer.BlockCopy(mValues, 0, result.mValues, 0, mValueCount * Constants.SizeOfDouble);
            Buffer.BlockCopy(lower.mValues, 0, result.mValues, mValueCount * Constants.SizeOfDouble, lower.mValueCount * Constants.SizeOfDouble);
        }

        private void StackSparseMatrices(SparseMatrix lower, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < lower.Rows; i++)
            {
                int begin;
                int end;
                lower.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i + Rows, lower.mColumnIndices[j]] = lower.mValues[j];
                }
            }
        }

        /// <summary>
        /// Concatenates this matrix with the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="right">The matrix to concatenate.</param>
        /// <param name="result">The combined matrix.</param>
        /// <exception cref="ArgumentNullException">If right is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Rows != right.Rows</strong>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not this.Rows x (this.Columns + right.Columns).</exception>
        public override void Append(Matrix right, Matrix result)
        {
            if (right == null)
            {
                throw new ArgumentNullException("right");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (right.Rows != Rows)
            {
                throw new NotConformableException("right", Resources.ParameterNotConformable);
            }

            if ((result.Rows != Rows) || (result.Columns != (Columns + right.Columns)))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = right as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                AppendSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    AppendSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < right.Rows; i++)
                    {
                        for (int j = 0; j < right.Columns; j++)
                        {
                            result[i, j + Columns] = right[i, j];
                        }
                    }
                }
            }
        }

        private void AppendSparseMatricesIntoNewSparseMatrix(SparseMatrix right, SparseMatrix result)
        {
            result.mValues = new double[mValueCount + right.mValueCount];
            result.mValueCount = mValueCount + right.mValueCount;
            result.mColumnIndices = new int[mValueCount + right.mValueCount];

            int rowStart = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                int begin;
                int end;

                FindRowStartAndEnd(i, out begin, out end);
                int length = 0;
                for (int j = begin; j < end; j++)
                {
                    result.mValues[rowStart + length] = mValues[j];
                    result.mColumnIndices[rowStart + length] = mColumnIndices[j];
                    length++;
                }

                right.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result.mValues[rowStart + length] = right.mValues[j];
                    result.mColumnIndices[rowStart + length] = right.mColumnIndices[j] + Columns;
                    length++;
                }

                result.mRowIndices[i] = rowStart;
                rowStart += length;
            }
        }

        private void AppendSparseMatrices(SparseMatrix right, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < right.Rows; i++)
            {
                int begin;
                int end;
                right.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, Columns + right.mColumnIndices[j]] = right.mValues[j];
                }
            }
        }

        /// <summary>
        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
        /// </summary>
        /// <param name="lower">The lower, right matrix.</param>
        /// <param name="result">The combined matrix</param>
        /// <exception cref="ArgumentNullException">If lower is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
        public override void DiagonalStack(Matrix lower, Matrix result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if ((result.Rows != (Rows + lower.Rows)) || (result.Columns != (Columns + lower.Columns)))
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = lower as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                DiagonalStackSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    DiagonalStackSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < lower.Rows; i++)
                    {
                        for (int j = 0; j < lower.Columns; j++)
                        {
                            result[i + Rows, j + Columns] = lower[i, j];
                        }
                    }
                }
            }
        }

        private void DiagonalStackSparseMatricesIntoNewSparseMatrix(SparseMatrix lower, SparseMatrix result)
        {
            result.mValues = new double[mValueCount + lower.mValueCount];
            result.mValueCount = mValueCount + lower.mValueCount;
            result.mColumnIndices = new int[mValueCount + lower.mValueCount];


            // Copy the row data
            Buffer.BlockCopy(mRowIndices, 0, result.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mRowIndices.Length; i++)
            {
                result.mRowIndices[i + mRowIndices.Length] = mValueCount + lower.mRowIndices[i];
            }

            // Copy the column data
            Buffer.BlockCopy(mColumnIndices, 0, result.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mValueCount; i++)
            {
                result.mColumnIndices[i + mValueCount] = Columns + lower.mColumnIndices[i];
            }

            // Copy the values
            Buffer.BlockCopy(mValues, 0, result.mValues, 0, mValueCount * Constants.SizeOfDouble);
            Buffer.BlockCopy(lower.mValues, 0, result.mValues, mValueCount * Constants.SizeOfDouble, lower.mValueCount * Constants.SizeOfDouble);
        }

        private void DiagonalStackSparseMatrices(SparseMatrix lower, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < lower.Rows; i++)
            {
                int begin;
                int end;
                lower.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i + Rows, Columns + lower.mColumnIndices[j]] = lower.mValues[j];
                }
            }
        }

        /// <summary>
        /// Returns a deep-copy clone of the <b>Matrix</b>.
        /// </summary>
        /// <returns>A deep-copy clone of the <b>Matrix</b>.</returns>
        public override Matrix Clone()
        {
            return new SparseMatrix(this);
        }
    }
}